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Abstract: Cosmological perturbations of sufficiently long wavelength admit a fluid dy¬ 
namic description. We consider modes with wavevectors below a scale k m for which the 
dynamics is only mildly non-linear. The leading effect of modes above that scale can be 
accounted for by effective non-equilibrium viscosity and pressure terms. For mildly non¬ 
linear scales, these mainly arise from momentum transport within the ideal and cold but 
inhomogeneous fluid, while momentum transport due to more microscopic degrees of free¬ 
dom is suppressed. As a consequence, concrete expressions with no free parameters, except 
the matching scale k m , can be derived from matching evolution equations to standard 
cosmological perturbation theory. Two-loop calculations of the matter power spectrum 
in the viscous theory lead to excellent agreement with iV-body simulations up to scales 
k = 0.2/i/Mpc. The convergence properties in the ultraviolet are better than for standard 
perturbation theory and the results are robust with respect to variations of the matching 
scale. 
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1 Introduction 

The evolution of cosmological perturbations and their growth at late times is a problem 
of fundamental significance in cosmology. The small amplitude of the primordial density 
perturbations over the isotropic and homogenous background makes linear perturbation 
theory a good description at early times. However, at late times the gravitational instability 
leads to growth of perturbations, inducing the emergence of large-scale structure and the 
formation of galaxies and galaxy clusters. As a result, linear perturbation theory ceases 
to be an accurate framework. Going beyond a linearized approach becomes mandatory for 
a precise comparison of theoretical predictions and astrophysical observations on length 
scales of 0(100Mpc) and below, where density and velocity perturbations are sufficiently 
large to render the evolution non-linear but still amenable to perturbative approaches. 
This range of scales is interesting because it displays the imprint of the oscillations of the 
primordial plasma on the matter density field, the so-called baryon acoustic oscillations 
(BAO). In this work, we show that viscous fluid dynamics provides a simple and accurate 
framework for extending the description of cosmological perturbations into this (mildly) 
non-linear regime. Compared to other formulations that we review below, our approach is 
novel in that it fixes the (effective) viscous coefficients and the pressure from a perturbative 
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calculation and without invoking any fit to non-perturbative information such as results 
from IV-body simulations. 

The main analytical approach that has been developed to go beyond linear theory aims 
at accounting for interactions between cosmological perturbations of different wavelength 
perturbatively in the amplitude of the primordial perturbations [1-4], see [5, 6] for reviews. 
This standard perturbation theory (SPT) is successful in producing the first non-linear 
corrections, but it is not a systematic approach for accurate calculations even for mildly 
non-linear scales. In particular, its convergence at the two-loop level becomes problematic 
at scales around 30-50 Mpc and the three-loop contribution becomes comparable to the 
linear spectrum at even larger length scales [7]. The physical origin of the apparent lack 
of convergence is that higher orders in the perturbative expansion become more and more 
sensitive to ultraviolet (UV) modes whose perturbative treatment is questionable. This 
has triggered the development of approaches beyond SPT aiming at a better control of 
perturbation theory at scales up to 10 Mpc [7-17]. Gaining a quantitative understanding 
of the BAO range motivates also the present work. 

One recent approach that aims at addressing this issue borrows ideas from the frame¬ 
work of effective field theory [18-26]. The suggestion is to treat the equations describing 
the large-scale perturbations as reflecting effective properties of dark matter that arise 
from integrating out the short-distance modes up to some coarse-graining scale. The im¬ 
plementation of this essentially Wilsonian point of view is complicated by the explicit time 
dependence of the process of growth of perturbations, and the resulting effective theory in 
general is non-local in time. Despite these difficulties, the formalism has been developed 
at the two-loop level and it compares well to the matter power spectrum of simulations 
down to scales of around 10 Mpc [25, 26]. The counterterms in this effective field theory 
of cosmological perturbations have the purpose of removing the dependence of loop cor¬ 
rections on the UV cutoff. They enter the effective energy-momentum tensor in the form 
of an effective pressure and viscosity coefficients, and their residual finite contribution to 
the spectrum is fixed by comparing with V-body simulations. 

Here, we follow a complementary approach to the growth of cosmological perturbations 
by starting from a viscous fluid dynamic formulation. This approach is motivated by the 
fact that essentially all physical systems lend themselves to a (not necessarily ideal) fluid 
dynamic description on sufficiently long time scales and for sufficiently long wavelengths 
[27]. Viscous fluid dynamics is an extension of ideal fluid dynamics that accounts for 
deviations from this idealization in a gradient expansion. As we shall discuss in detail 
in the following, the viscous fluid dynamic approach to large scale structure has some 
commonalities with the effective field theory approach discussed above. In the present 
work we focus on the physical origin of the non-ideal behavior on large scales. For this, it 
is conceptually important to realize that applications of viscous fluid dynamics allow for 
viscosities of two qualitatively different origins: 

• Fundamental (microscopic) viscosity: This viscosity is due to the diffusive transport 
of momentum by the microscopic constituents of the fluid, that is, by particles or 
radiation. It is a fundamental property of the matter under consideration that can 
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be traced back to the fundamental interactions between the fluid constituents. In 
cases in which the quantum field theory of the matter is known, these fundamental 
viscosities are calculable from the Lagrangian via the Green-Kubo formula [28, 29]. 

• Effective viscosity: If one limits a fluid dynamic description to large wavelengths, 
i.e. wavenumbers k below a certain matching scale k m , then dissipation of these fluid 
perturbations can be mediated by fluid dynamic degrees of freedom on scales k > k m . 
The use of effective viscosities parametrizes then how momentum carried by the long 
wavelength modes dissipates to shorter wavelengths (k > k m ) that are fluid dynamic 
in nature but that are not followed in the large wavelength fluid dynamic description. 1 
The physical mechanism underlying this is the interaction between long and short 
wavelengths modes. Examples are ubiquitous and include e.g. the description of 
approximately laminar flow on large scales in the presence of short scale eddies (so- 
called eddy viscosity) [30]. The introduction of effective viscosities is a method of 
choice, in particular if one aims at limiting an application of fluid dynamics to scales 
that evolve with only mildly non-linear effects (where semi-analytical methods apply), 
while smaller length scales show a strong non-linear (even turbulent) behavior with 
certain degrees of freedom even outside of the fluid approximation. In contrast to 
fundamental viscosity, effective viscosity is not a fundamental property of the matter, 
but it depends on the spectrum of excitations k > k m with which the long-wavelength 
modes can interact, and on the matching scale k m up to which these long-wavelength 
modes are evolved explicitly in the fluid dynamic description. 2 

We emphasize that the viscous fluid dynamic formalism introduced in section 2 allows 
for exploring observable consequences of dark matter with either fundamental or effective 
viscosities and pressure. The fundamental (microscopic) material properties of the cosmo¬ 
logical fluid at late times (dominated by dark matter) may have observational consequences 
on galaxy cluster scales [31-36], as well as on very large cosmological scales [37-44], Still it 
is well conceivable that microscopic viscosity and pressure are negligible for the evolution 
of cosmological perturbations up to the BAO range. In particular, they may be much 
smaller than the effective viscosity and pressure of a fluid dynamic description limited to 
k < k m ~ 0(lh/Mpc). An analogous statement holds for the notion of pressure: dark 
matter may be essentially pressureless on a fundamental level, while an effective pressure 
can arise if fluid dynamics is limited to k < k rn . 

In the present work, we explore specifically the possibility that effective viscosities and 
pressure dominate for the description of scales k < k m ~ O (1 h /Mpc). For this case, we 
argue in section 3 that realistic numerical proxies for the effective viscosities and pressure 
can be obtained by matching viscous fluid dynamics to results from standard perturbation 

1 The fluid description is expected to break down on very short length scales. The corresponding modes, 
in principle, can affect the long scale dynamics. However, as discussed in section 3, there is evidence that 
this effect is small on BAO scales. 

2 Microscopic viscosity could, however, play a role at short wavelength. There could also be an interplay 
between the effective and the microscopic viscosity, in the sense that momentum is transported to the short 
wavelength fluid modes by the former, from where it is then dissipated by the latter. 
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theory. In the same section 3, we also discuss issues about the stability of cosmological 
standard perturbation theory and to what extent our matching procedure is independent 
of those. In section 4, we then explore numerical solutions of viscous fluid dynamics based 
on this matching procedure, and we show that this gives rise to a satisfactory description of 
large scale structure and its red-shift dependence in the BAO range up to scale ~ 0.2 h/ Mpc. 
We finally summarize our main findings in the conclusions. 

2 Viscous fluid dynamics and cosmological perturbations 

The evolution of cosmological perturbations is governed by the Einstein equations 

Gp,u = 8 ttGnT . ( 2 - 1 ) 

These equations imply the conservation of the energy-momentum tensor 

V U T^ U = 0 , ( 2 . 2 ) 

which for the case of an ideal fluid without any conserved charges completely defines its 
dynamics. Here, we seek a more realistic (generic) description of the cosmological medium 
by considering viscous-fluid dynamics as an effective theory of cosmological perturbations 
of sufficiently long wavelengths. This approach goes beyond the perfect fluid approximation 
by accounting for deviations from local equilibrium up to first order in gradients, and it 
allows for non-vanishing shear and bulk viscosity. Imperfect fluids have been considered 
in the past for different aspects of cosmology at large scales (an incomplete list of relevant 
publications includes [37, 38, 41, 43, 45-49]). While the main effort in these works has 
been the implications for dark energy or linear behaviour of perturbations, we will be more 
interested in the growth of structure in the mildly non-linear evolution of dark-matter 
perturbations. 

The energy-momentum tensor for a relativistic viscous fluid is of the form 

= pu»u v + {p + ir b )A^ + 7 r^. (2.3) 

Here p is the energy density and p the pressure in the fluid rest frame, 717 , the bulk viscous 
pressure, and the shear-viscous tensor, satisfying: = 7 r /J /t = 0. The matrix A^ 

projects to the subspace orthogonal to the fluid velocity: A^ u = g^ u + u^u v . The bulk 
viscous pressure is 

7T b = -CVX, ( 2 - 4 ) 

where ( is the bulk viscosity, while the shear viscous tensor is given by 

7f lw = -2 rja^ = —2rj Q (A'"V^ + A ra V Q u") - ^A^(V«a a )^ , (2.5) 

with shear viscosity p. We work in the first-order formalism of relativistic fluid dynamics 
and drop terms of second order in gradients. In the non-relativistic limit this approxi¬ 
mation leads to the standard Navier-Stokes theory. This should be sufficiently accurate 
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for cosmological perturbations well inside the horizon and that relaxation times and other 
second-order terms would only lead to minor quantitative modifications, although they are 
in principle needed for a viable causal structure and for linear stability [50, 51]. 

The viscous fluid dynamic equations (2.2) take then the explicit form 

+ (p + - C (V M u M ) 2 - 2 tja^a^ = 0 , (2.6) 

(p + P + 7 T b )u^^u a + A^V M (;p + v r h ) + = 0 • (2.7) 

The Einstein equations (2.1) connect perturbations in the matter fields to metric per¬ 
turbations. The latter enter (2.6), (2.7) via the covariant derivative V^, the projector 
A / "' / = g^ u + u^u 11 , and the form of the four-velocity u^ = dx^/y/—ds 2 . For the metric, we 
consider an ansatz of the form 


ds 2 = a 2 (r) [—(1 + 2'h(r, x)) dr 2 + (1 — 2<h(r, x)) dxdx] , (2.8) 


that accounts for the dominant scalar metric perturbations only. These are parametrized by 
the Newtonian potentials $ and 'll. From the linearized spatial components of the Einstein 
equations, one can check that the difference — T is proportional to shear viscosity, 
see eqs. (2.19), (2.20) below. Keeping two Newtonian potentials is thus necessary for 
the description of a fluid with non-vanishing shear viscosity. Since the flow field u 11 is 
normalized, = —1, it can be expressed through the coordinate velocity v l = dx l /dr 

[i,j = 1, 2,3) and the potentials <h and \k: 


U a^J 1 + 2T- (l-2$)u 2 ^ 1 ’ t; ' ) 

In the next subsection we turn to an analysis of the fluid dynamic equations to linear order 
in perturbations. Within this framework, we discuss in subsection 2.2 how the growth of 
density perturbations is modified by pressure and viscosity. In subsection 2.3, we extend 
our analysis to the dominant non-linear contributions at the scales of interest. 


2.1 Linear analysis 

We parameterize the density and pressure in terms of spatially homogeneous and isotropic 
background fields p m (r), p m (r) and small spatially varying perturbations, p(r, x) = p m (r)+ 
bp(r, x) and p(r, x) = p m {r) +5p(r, x), with 5p, bp <C p m ■ We consider a flat Universe with 
cosmological constant A and one species of viscous matter. With the help of the definitions 

n = a , H = -n, (2.io) 

a a 

the equations for the background fields are 

7U = —2 = ^Gat(A + p m )a 2 , (2-11) 

a z 3 

Pm + 3 (p m + Pm)7d = 0 . ( 2 - 12 ) 

We consider perturbations in the variables bp, bp, v l , d>, T up to linear order, and we focus 
on scalar perturbations only. To expand the fluid dynamic equations in dimensionless 
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measures of these perturbations, we use the normalized density contrast 5 and the velocity 
divergence 9 , 


8 = 


5p_ 

Pm 


9 = W. 


(2.13) 


Also, 5p is not an independent perturbation, but is related to density perturbations by the 
equation of state. Specializing to a simple equation of state, we write 


Pm 

5 

Pm 



Cad 


-Pm _ 

Pm 


W 

3(1 + w) 


n, 


(2.14) 


where c s is the velocity of sound. Finally, we introduce shorthands for the normalized 
so-called kinematic viscosities 


v = 


( Pm+Pm)a (1 + w)p m a ’ 


"C = 


C 


(1 + w)p m a 


(2.15) 


For notational economy, the dependence of the evolution equations on bulk viscosity 
is not made explicit in our discussion. However, it can be restored by remembering that 
bulk viscosity £ enters the linearized fluid dynamic equations in only two ways. First, it 
contributes to the effective pressure of the system, where its contribution can be restored by 
replacing p m —»• p m — 3 ( H in our expressions. 3 Second, bulk viscosity arises in terms that 
describe the dissipative attenuation of perturbations. For scalar perturbations it appears 
in the combination and can be made explicit by replacing |zz —y With 

these provisos, evolution equations for cosmological perturbations in bulk and shear viscous 
matter can be written in Fourier space as 


4 = -(1 + w) (o k - 3<b k ) - 3 U{c 2 s - w)5 k , (2.16) 

4 = ~U9 k + kH k + + 3 c 2 ad H9 k - \k 2 u9 k . (2.17) 

Here, we have introduced the Fourier-transformed fields 

8 ( t , x) = J d 3 k e ,kx 4(r), 6 ( t , x) = J d 3 k e lkx 0 k (r). (2-18) 


The equations of motion (2.16), (2.17) can be closed by the Poisson-like equations for the 
Newtonian potentials, 

k 2 $ k = -\si m U 2 ^3(1 + ™)^^ + 4) , (2.19) 

k 2 ^ k = ——^3(1 + w )~^2^ + 4 + 4zz0 k ^ • (2.20) 

3 If one allows for a density dependence of (,, this bulk viscous pressure modifies not only the background 
equation (2.12), but also the derivatives in (2.14). We note that, in principle, p m — 3 ( H can be negative, 
see e. g. ref. [41] and references therein. Our formalism can be adapted to this case, but we do not discuss 
this point further. 
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For the case of a pressureless ideal fluid, w = c 2 = c 2 d = 0 and v = 0, there is 
only one Newtonian potential, and the equations (2.16), (2.17) simplify. To understand 
these simplifications and how they extend to a non-ideal fluid, we discuss now the scales 
in the problem. The time-scale of the evolution is set by the scale 77, and we consider 
subhorizon perturbations with wavenumbers k 77. A time derivative is thus equivalent 
to a multiplicative factor of 77, and it is consistent to assume that the Fourier-transformed 
fields 5k and 8 k /7~L are of comparable magnitude. Eqs. (2.19), (2.20) indicate that 4>k and 
\kk are suppressed by 77 2 /k 2 relative to 5'k and 6 k /'H. The above equations then read 


4 = -8 k + 0 (7i 2 /k 2 ) , 

(ideal fluid) 

(2.21) 

e k = -w k + kH k + o in 2 /k 2 ), 

(ideal fluid) 

(2.22) 

k 2 $ k = A; 2 ^ = ~ln m n 2 8 k + o (' u 2 /k 2 ). 

(ideal fluid) 

(2.23) 


To study within this linear analysis the late time evolution of the matter distribution, we 
shall allow 5k and 0 k /'H to take values up to order 1. On the other hand, it is apparent 
that a perturbative analysis cannot be valid for 5k, 9 k /7-L S> 1. 

The presence of non-zero pressure and viscosity can disrupt the hierarchy of scales that 
we assumed. In particular, for c 2 ,w;,z/77 of order 1, the terms ~ k 2 8 k and ~ fc 2 $k would 
dominate in eq. (2.17). They would result in an evolution with a characteristic times-scale 
much shorter than 1 /77, related to the appearance of density waves or exponential damping 
of fluctuations through viscosity. In principle, such strong attenuations are conceivable. 
Here, however, we aim at maintaining an effective description with 1/77 as the only relevant 
timescale. This corresponds to the assumption that w, c 2 , c 2 d , ul~L ~ 77 2 /A/( l , where k m is 
the largest wavenumber (corresponding to the shortest length scale) for which this effective 
description is applicable. To make this point explicit, we introduce the parametrization 


c 


2 

s 


vH 


a & 


77 ^ 

U2 ’ 

form 


3 77 2 

1 au W 


(2.24) 

(2.25) 


and we assume that a s , a v are at most of order 1. We shall confirm that this assumption 
is consistent with the values obtained by matching viscous fluid dynamics to results from 
standard perturbation theory, cf. sec. 3. We shall also provide numerical evidence that 
values of this order provide a good description of the power spectrum in the BAO range. 
The linearized equations then become 


4 = ~0 k + O (U 2 /k 2 ) , (2.26) 

8 k = ~H9 k + kH k + a s ^k 2 5 k - k 2 aj^ 6 ^ + O (U 2 /k 2 ) , (2.27) 

Km K m rL 

k 2 $ k = fe 2 T k = -^ m -H 2 5 k + O ( U 2 /k 2 ) , (2.28) 

where we have kept again terms up to O ('H 2 /J‘: 2 ) only. We note that up to this order, 
effects from finite pressure other than parametrized by the velocity of sound c 2 do not enter 
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since they are not enhanced by factors k 2 in equations (2.16) and (2.17). In the following, 
we sometimes refer to the sound velocity term as effective pressure. 

We are now in a position to discuss to what extent effective viscosity and pressure 
terms can delay or inhibit the growth of structure. To this end, we combine the linearized 
equations (2.26), (2.27) and (2.28) into a single second-order differential equation for the 
density contrast <5k 

<5k + i^T~L + 

Changing the evolution variable to rj = In a, this can be written as 

— ~ ° is W~ S ) ^' k = 0 > (2.30) 

where the prime denotes a derivative with respect to 77. 4 


k z u ) <5 k - [ 'Jklm'H 2, - elk 2 ) (i'k = 0. 


(2.29) 


2.2 The linear growth factor in the presence of viscosity and pressure 

The solution of (2.30) defines the linear growth factor. In the absence of effective pres¬ 
sure and shear viscosity, this growth factor is independent of the wavelength 2n/k of the 
perturbation. For instance, for the simple case of ACDM cosmology, and neglecting the 
radiative components, we have %'/T-L = 1 — 30 m /2, and, for a v = a s = 0, the amplitude 
of the solution takes the well-known form 
/H(a) f a da' 


D L (a ) = A- 


(for ACDM without radiative component) (2.31) 


u Jo 

with A a normalization constant. 

For non-zero effective pressure and shear viscosity, a u ,a s / 0, eq. (2.30) shows that 
the growth of cosmological perturbations depends on their wave vector k [43, 45, 47, 49]. 
Moreover, since the effective pressure and viscosity will depend in general on the matter 
density (and possibly on other characteristics), a v = a u (p m ), a s = a s (pm), we expect a 
non-trivial scale dependence a u = a u (rj), a s = en s {p) that reflects the (effective) properties 
of dark matter. 

The qualitative behavior of the solutions of eq. (2.30) can be understood by the ana¬ 
lytical study of simplified set-ups. To this end, we consider an Einstein-de Sitter Universe 
with Q m = 1, for which eq. (2.30) reduces to 


with 


5Z 


A u = 


k /^Asd k — 0, 

(2.32) 

k 2 

+ Of,p-, 

(2.33) 

(3 k 2 \ 

U as ki)- 

(2.34) 


4 Note that we follow the conventional notation p = ln(a), which should not be confused with the shear 
viscosity introduced earlier. We parameterize viscosity by the coefficient a v defined in (2.25) and (2.15) in 
the following. 
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For constant values of a v and a s , we obtain a solution of the form 



where c\ and C 2 are integration constants. In the long wavelength limit k/k m —> 0, where 
effective pressure and viscosity become unimportant, one recovers from this solution the 
standard growing mode ~ exp(? 7 ) and the decaying mode ~ exp(—3r//2). Aside from this 
consistency check, the wavelength dependence of (2.35) has several interesting features: 

1. The case of finite effective pressure without effective viscosity (a s > 0, a v = 0/: 

For increasing wavenumber k, A s decreases, and this reduces the exponential growth 
to ~ exp (dp) with d < 1. For k 2 /k^ > 3/(2 a u ), A s becomes negative, and for the 
slightly higher wavenumber k 2 /k^ > 25/(16 a v ), the argument in the square root 
of (2.35) becomes negative and the solutions become oscillatory. This is what one 
expects from pressure: on scales on which pressure becomes important, it counter¬ 
acts the growth of the density contrast due to gravitational collapse; the system 
bounces back and starts oscillating. Wavelengths for which this happens will not 
form structure, since the real part of both modes in (2.35) is decaying. 

2. The case of finite effective viscosity without effective pressure (a s = 0, a v > 0/: 

For the mode ~ exp {dr\) with d = \ A v + \JJff + , the exponential factor d 

takes the value of the standard growing mode d = 1 for vanishing viscosity and it 
approaches zero as d 3/(1 + 2a v k 2 /kff) > 0 in the limit of very large wavenumber 
or viscosity. This mode is always growing, d > 0. Thus, in marked contrast to the 
case of finite pressure, finite viscosity does not inhibit structure formation at any 
wavelength, but it slows down the formation of structure on small scales. 

3. The case a s > 0, a u > 0: 

For perturbations with k 2 /k 2 n > 3/(2cc s ), A s becomes negative and both modes 
of the solution (2.35) are decaying. Thus, in a viscous medium as well, pressure 
inhibits structure formation on small length scales. However, for perturbations that 
are sufficiently far in the UV, the square root y/A 2 + A s is now real, since the factor 
a u k 2 /k^ enters in quadrature: the oscillations that pure pressure induces are now 
overdamped by viscous effects. 

Analytical solutions of (2.32) can be given even for some time-dependent parameters 
a u , a s . Here, we consider 


a u = f3 v exp (nr]) 


a s = fi s exp (up) 


(2.36) 


for which the solution to (2.35) can be written in terms of the hypergeometric function 

iFi(a, b, z): 



(2.37) 
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Of particular interest for the following is the case k = 2, since we shall motivate in section 3 
an ansatz with a time-dependence that reduces to eq. (2.36) for this value of n. Fig. 1 shows 
the solution (2.37) for n = 2 and a set of wavenumbers k. One sees that in comparison 
to the pressureless non-viscous case, the linear growth factor is reduced due to finite shear 
viscosity and pressure. This reduction is very significant for large wavenumbers (small 
scales), while for the parameters chosen in fig. 1 scales larger than ~ lOMpc are unaffected 
by the effective pressure and viscosity. 



Figure 1. The time- and wavenumber-dependence of the linear growth factor D v i sc (k,z) (given 
by (2.37) for ci = 1, C 2 = 0) for an Einstein-de Sitter cosmology in which dark matter is endowed 
with an effective pressure and viscosity, normalized by the corresponding growth factor Dl(z) for 
pressureless non-dissipative dark matter (being identical to the scale-factor a = e v in Einstein-de 
Sitter). Shown is the time-dependence for various wavenumbers k assuming ol s ,ol v ~ exp(2 77 ), 
with k m = 0.6/i/Mpc and j3 s = (3 U = 0.67 (this particular choice will be motivated later in 
section 3) . The result also applies to ACDM cosmologies to a very good accuracy when substituting 
77 —> fj = In (Dl(z)), and a* —> on, see eq. (2.57). 


The main features of the solution (2.37) for k = 2 are similar to those of the solution 
(2.35) for time-independent a s , a v . In particular, for the case of finite pressure ({3 S > 0) 
but vanishing viscosity (f} v = 0), the solution of eqs. (2.30), (2.36) simplifies to 


4 (??) 


ci exp(— 77 /4) J _ 5/ 4 



+ ci exp(— 77 / 4 ) J 5/4 



.(2.38) 


Expanding in the argument of the Bessel functions for small /3 S (which at fixed time 77 
is a fortiori an expansion for sufficiently long wavelength 2n/k), one recovers the stan¬ 
dard growing and decaying mode for /3 S —> 0. On the other hand, it is apparent from 
(2.38) that for sufficiently short wavelength 2n/k, both modes are decaying. Thus, pres- 
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sure counteracts structure formation on small wavelengths, in close similarity to the time- 
independent case discussed above. For the case of finite viscosity (f3 u > 0) but vanish¬ 
ing pressure (j3 s = 0), one can expand (2.37) for large viscosity (or equivalently large 
k). For the growing mode one finds in this limit <5^ -» c\\j2f/3 U [r(|)/r(|)] (k m /k). 
Thus the growth is slowed down and asymptotically the linear solution saturates at a 
finite, scale-dependent value at late times. This result can be generalized to k / 2, 
<5k —> ci [T(l + 5/2a)/r(l + 3/2a)] (k m /k) 2 ^ K . We conclude that in marked con¬ 

trast to the physical mechanism induced by pressure, finite viscosity does not inhibit struc¬ 
ture formation, but it slows it down. 

In the opposite limit of small viscosity/pressure, or equivalently large wave length 
2ir/k, the leading corrections are for n = 2 given by 

4\ 


M 7 ?) -* °i exp(77) 1 - 


OL v + OL s 

9 


k 

TT 

h ID 


2 2>af y + 4 a u a s + a 2 


234 


k_ 

km. 


(2.39) 


+ C 2 exp(—3?y/2) 1 - 


2Of e 


2 3a 2 + 4 a u a s — 4a 2 

+ 24 



As expected, the leading correction scales as k 2 , and depends only on the sum a v + a s 
for the growing mode. This degeneracy is broken by the k 4 corrections, as well as by the 
decaying mode. 

In summary, while the analytical structure of (2.37) is more complex and thus more 
difficult to discuss than that of the solution (2.35) for time-independent viscosity and 
pressure, our short discussion of (2.37) points to features of the role of viscosity and pressure 
that are model-independent. 


2.3 The dominant non-linear terms 

The complete covariant fluid equations (2.6) and (2.7), supplemented by equations for the 
metric (Einstein’s equations) contain many non-linear terms. To identify amongst them 
the dominant ones for the problem at hand, we take recourse to the hierarchy of scales 
suggested by the linear analysis in section 2.1. More specifically, we make the following 
assumptions: 

• We rely on perturbation theory. This means that, during manipulations of the evo¬ 
lution equations, we assume that the fields 5 = 5p/p m , O/kL and the potentials are 
small. However, we allow also for values of order 1 at the limit of validity of our 
analysis. 

• With respect to H/k, we treat 5 and 6/T-L as quantities of order 1, v as a quantity of 
order %/k and T, <F as quantities of order H 2 /k 2 . 

• We assume that a time derivative is equivalent to a factor of 77, while a spatial 
derivative to a factor of k. 

• We assume that the couplings w , c 2 , c 2 d , vH are of order 'H 2 /kf n , where k m is the 
largest wavenumber (corresponding to the shortest length scale) for which this effec¬ 
tive description is applicable. 
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Within this hierarchy, the linear evolution equations given in eqs. (2.26), (2.27), (2.28) 
up to order H 2 /k 2 receive the following dominant non-linear corrections 


4 + 4 + / d 3 p d 3 q <5 (3) (k - p - q) cq(p,q) <5 p 6> q = 0 


4 + ( H + ^uk z ) 6 k + ( ~-Sl m U z - c^k 2 ) 4 , 


q) s p 6» q = 0, 

(2.40) 

14, 

(2.41) 


+ / d 3 pd 3 q5 (3) (k-p-q) (/?i(p,q)44 + /3 2 (p,q)6 , p 6» q + /? 3 (p,q)4 5 6 'q) = 0 ’ 


where ' J 

ai(p,q) 

/?i(p,q) 
P 2(p,q) 

/?3(p, q) 


(p + q)q 

9 ’ 

<r 

(2.42) 

c 2 (p + q)q, 

(2.43) 

(p + q) 2 pq 

2 p 2 q 2 

(2.44) 

4 / N 
-~Kp + q)q- 

(2.45) 


A particularly useful, compact formulation can be obtained by rewriting these equations 
for the doublet 


^</>i(??,k) 

\ 02(?7, k) 


/ 44) \ 

V n ) 


(2.46) 


and the evolution parameter rj = In a(r). One obtains 


<4<4(k) = -n ab (k,r/)(j) b (k) + j d 3 pd 3 q6^ 3 \k - p - q) 7 a & c (p, q, rj) <f> b (p) <j> c {q), (2-47) 


where 


^(k,r?) = 


o “1“ ^ 


s 1.2 


1 + 


- 1 
- + ex 


k*_ 

V U2 


(2.48) 


where the prime denotes a derivative with respect to 77 . The evolution equation (2.47) is 
complete up to order H 2 /k 2 The non-zero elements of 7 a bc are 


7i2i(q,p,??) = 7ii2(p,q,7) 
7211 (p,q,??) 
7222(P,q,?i) 
72i2(p,q,??) = 722i(q,p,7) 


«i(p^q) = (p + q)q 
2 2 q 2 ’ 

/3i(p,q) _ (p + q)q 

H 2 as kl 


/4(p,q) 


(p + q) 2 p • q 

2 p 2 q 2 


A(p,q) _ (p + q)q 
2 n av 2 k 2 m 


(2.49) 

(2.50) 

(2.51) 

(2.52) 


5 Depending on whether we define dimensionless thermodynamic variables by normalizing with respect 
to pm , as done in eqs. (2.13), (2.14) and (2.15), or with respect to p m + 5p, the structure of the couplings 

(2.43), (2.45) will change. We shall demonstrate in section 4 that this ambiguity is numerically unimportant. 
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In an Einstein-de Sitter Universe, the linear growth factor is Dl = a and the time 
variable 77 in (2.47) can thus be interpreted as In Dl- In general, this is not the case, but it 
turns out to be convenient to rewrite the equations of motion with a modified time variable 


fj = In D l , 


For the suitably redefined doublet 

^i(r),k) 

k) 

the evolution equations take the form 


dp _ 1 dP L 
dp Dl dp 

( 4(t) > 

_^k(r) 

V Uf / 


(2.53) 


(2.54) 


<9^ a (k) = - ^^(k, 77 ) + <50(k, 77 )^ 4> b (k) + j d 3 p d 3 q <5(k — p — q) 7 abc(p, q, v) 4>b(p) 4>c(<l), 

(2.55) 

with 

n L {k,fj) =( 3 J m 3 n m 1 1 ) : <50(k, fj) = ( fc °2 ° k 2 J , (2.56) 

\~ 2 JT 573 -- 1 J a ^J 

and 

Of, = jO v , a s = J 2 ’ av = ~f~- (2.57) 

The vertices 7 a bc are identical to 7 a bc except for the replacement a* —> dj. The advan¬ 
tage of this formulation is that, even for the ACDM cosmology, we have to a very good 
approximation 

li L (k,i))~f 3 jV (2.58) 

This is the form of the leading order expression in an Einstein-de Sitter Universe. The 
transformation to the time variable p will thus allow us in the following to adapt perturba¬ 
tive results derived for an Einstein-de Sitter Universe to a discussion of ACDM cosmology. 


3 Matching effective pressure and viscosity to results from standard per¬ 
turbation theory 


To describe large scale structure with the viscous fluid evolution equations (2.55), (2.56), 
we need to specify the functions a s (rj), a u (fj) and the scale k m entering these equations. 
In the present section, we do this by matching the fluid dynamic evolution to results for 
the effective propagator [9] 


/ 5(j>a(k,f)) \ 

W(k',7')/ ’ 


Gab(k, fj, rj')5^ (k — k') 


(3.1) 








calculated from standard perturbation theory beyond linear order. This propagator defines 
the response of the field at late times r) to a perturbation at earlier times ff. Before pre¬ 
senting technical details of this matching, we discuss first the motivation for this procedure: 

As mentioned already in the introduction, viscous fluid dynamics accounts for a wide 
variety of physically realized systems. The viscosities used in these descriptions are not 
necessarily the fundamental ones, by which we mean that they do not reflect diffusive 
transport mediated by the fundamental constituents of the liquid. Rather, the dominant 
source of diffusion is often mediated by fluid dynamic degrees of freedom. This is the case in 
particular if one limits a fluid dynamic description to long wavelengths, while perturbations 
of short wavelengths (that are in principle fluid dynamical and whose evolution can be 
highly non-linear) are not followed explicitly. These UV modes modify the evolution of the 
IR modes. This modification can be accounted for by introducing scale-dependent effective 
viscosities and pressure in the fluid dynamic description of the long-wavelength modes. 
Examples for this procedure are ubiquitous. They include, for instance, the effective eddy 
viscosity in cases where laminar flow on large scales is modified due to the presence of 
small-scale turbulent eddies in the system [30]. There are also examples in which the non¬ 
linear UV modes are not vortical degrees of freedom (e.g. wave turbulence). Here, we seek 
a viscous fluid description of the cosmological fluid on sufficiently long wavelengths without 
following the evolution on short scales. To this end, we aim at choosing a s (fj), ck v (rj) as 
effective viscosity and pressure in the sense described above. In contrast to fundamental 
properties of matter, these effective properties depend on the scale up to which the effective 
fluid dynamic description is applied. They also depend on the spectrum of perturbations 
above this scale. The matching proposed below will make these features explicit. 

The concept of an effective theory of the low-A; modes of the cosmological fluid [18— 
22, 24, 25] is closely related to the viscous-fluid dynamic description formulated here. This 
low -k effective theory contains additional terms that are not present in the fundamental 
description of dark matter. From a Wilsonian point of view, these terms account for the 
effect of the large-/c fluctuations that are integrated out in the coarse-grained theory. In 
close similarity to the effective material properties in the viscous-fluid dynamical descrip¬ 
tion, the couplings in the effective theory depend on the spectrum of UV perturbations 
and on the cutoff scale up to which the effect of these perturbations is coarse-grained. 

In principle, the magnitude and time-dependence of the effective pressure and viscosity 
can depend on the dynamics of UV modes in a complicated way. Therefore, in the approach 
followed in Refs. [19, 22, 25] these quantities were essentially treated as free parameters (or, 
more precisely, free functions of time) that have to be determined by a fit to V-body data. 
On the other side, using IV-body simulations it was found that the impact of the strongly 
non-linear physics at small scales on the power spectrum in the BAO range is relatively 
weak [55, 56, 58]. Motivated by this, we explore in the following the possibility that for a 
suitable choice of the coarse graining scale, the dominant contribution to effective viscosity 
and pressure is generated by modes for which a perturbative treatment within the fluid 
picture is still applicable. This set-up is formulated in the next subsection 3.1, where we 
match d s (f/), a u (fj) to results from one-loop perturbation theory. In subsection 3.2 we 
then discuss related issues of perturbative stability of standard perturbation theory and 
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comment on the possible impact of strongly non-linear dynamics at very short scales. 


3.1 A simple matching procedure to determine effective viscosity and pressure 


The calculation of the effective propagator (3.1) in standard perturbation theory provides 
a particularly simple example of how effects of UV modes modify the propagation of long- 
wavelength cosmological perturbations. Assuming that dark matter is an ideal pressureless 
fluid on the fundamental scale, and specializing to the limit of large separation between 
fj and fj' (so that only the fastest growing mode needs to be kept), standard perturbation 
theory yields for small k and in an Einstein-de Sitter Universe [9] 


G ab ( k, fj, fj') = g ab (fj - fj') + G^ b (k, fj, fj’) = g ab {fj - fj') - k 2 e v v 'crj{fj) 



(3.2) 


Although these results have been derived for the Einstein-de Sitter Universe, they can also 
be used for ACDM due to the held redefinition (2.54) for the approximation in eq. (2.58). 
Here, the linear propagator is 


9ab(fj - fj') 



e -3(»W)/2 

5 




(3.3) 


where the decaying mode in the second term becomes unimportant in the limit of large fj—fj' 
for which (3.2) is written. To lowest order, cosmological perturbations propagate linearly. 
To one-loop, however, the propagator (3.2) depends on the dimensionful parameter 

A POO 

<Z = Yj 0 dqP L (q), (3.4) 

where P L (k) stands for the density-density spectrum. The redshift dependence of a 2 (rj) 
in (3.2) arises from the redshift dependence of P L (q). Equations (3.2), (3.4) show that in 
standard perturbation theory at one-loop, a cosmological perturbation measured at time 
fj on scale k will be influenced by the spectrum P L {q) at all scales q. 

We seek a viscous fluid dynamic description of the cosmological fluid for scales k < k m 
only. Here k m is a matching scale scale chosen such that strongly non-linear perturbations 
have wavenumber larger than k m . Perturbations with k > k m will not be evolved explicitly, 
but their effect on the long wavelength modes will be parametrized in terms of the effective 
viscosity and pressure, a s (fj), a u (fj). Since the fluid dynamic equations of motion will 
account explicitly for the interactions of modes k < k m with modes q < k m , we must not 
absorb such interactions in the definition of a s (fj), a y (fj). These considerations prompt 
us to match the expression for the propagator in viscous fluid dynamics to the standard 
one-loop result (3.2), with the replacement 6 

POO 

°d —> o-rf,match = y / dqP L (q) . (3.5) 

J km 

Another formulation that used the propagator to ameliorate the non-linear predictions of SPT is renor¬ 
malized perturbation theory [8]. This proposal yields a description to resum the infrared modes of the 
theory and can reproduce the non-linear evolution of the BAO peak [16]. Notice that our method makes 
use of the opposite regime of UV scales, whose influence we treat in terms of the effective viscosities. 
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For simplicity, we use the same matching scale k m for all redshifts z. As a consequence, 
the independence of Od )ina t c h results only from the linear growth factor. This allows us to 
parametrize <7 d , match as 


2 _ P(km) Dl{v) 

° d ,match \ I) r\ 7,2 


(3.6) 


An estimate for /3(k m ) can be obtained by computing cr dmatch for a spectrum that scales 
exactly as k 3 . We find 


P(k m ) = —kl l P L {k m ). 


(3.7) 


Evaluating fi{k m ) from this expression results in values near 1 for k m in the range 0.4 — 
1 h/Mpc. For our numerical analysis we shall compute f5{k m ) from the exact linear spec¬ 
trum. 

Our matching procedure implies that we should regard (3.2) as the solution of an 
equation of motion for the long-time and long-distance behavior that can be written as 


( S ac dfj + H ac (k, fj)^ G cb ( k, fj, fj') = 0. (3.8) 

We are interested in the limit of this equation for small k. For k —> 0 the effective 
propagator approaches the conventional linear propagator and H(k, fj) —> Q L (fj). The 
leading correction is determined by substituting the one-loop propagator (3.2) in eq. (3.8) 
and solving for O(k, fj) to order k 2 . We find 


n(k, fj) 



^ + k 2 a 2 d {fj) 



(3.9) 


where we have used that &d(fj) = exp (77 — fj') 0 ^( 77 ') on the growing mode. We can now 
compare this expression to the form of the linear ff(k, fj) in the effective theory, as given 
by eq. (2.56). Apart from numerically small terms in the first row of the second matrix on 
the right hand side of eq. (3.9) that amount to ~ 1% corrections, we can match the two 
expressions if we take 

® s j 2 gg &d ,match vl) yQ P Dl\Z) , 

Uv = y = ^*4 °d, match (V) = f j D l(. z ) ■ (3-fO) 

Eqs.(3.10) yield separate values for a s and a u . However, their determination is based 
on the one-loop expression (3.2) for which only contributions from the growing mode have 
been kept. Beyond this approximation, the values in (3.9) are expected to change. To 
make ambiguities in this matching apparent, one may therefore project the fluid dynamic 
expression for SQ(k,fj) on the growing mode, 


hH(k, fj) 



~( 0 0 A 

5 l 3(d<j + djy)p- 2(a s + <5^)p- ) 

\ *"m / 


(3.11) 


and compare it to the second term on the right hand side of (3.9). This shows that the 
sum a s + ot u has greater significance than the separate values in (3.10). This sum matches 
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Figure 2. Linear matter power spectrum for a ACDM model with erg = 0.79, fi m = 0.26, fib = 
0.044, h = 0.72, n s = 0.96, and the one- and two-loop corrections in standard perturbation theory at 
z = 0. The grey shaded lines show analytical approximations valid for large or small k, respectively, 
see [52] for details. 

exactly the lower row of (3.9). The upper row of (3.9) cannot be matched in this procedure, 
but in comparison to the lower row it is only a 1% correction. 

3.2 The reliability of the perturbative spectrum and the effective viscosity 
and pressure 

In the previous subsection we fixed the effective viscosity and pressure of the fluid dynamic 
formulation of large scale structure through eqs. (3.5), (3.10). Here we would like to 
understand for which range of matching scales k m one expects this mapping to be accurate, 
and furthermore discuss various sources of uncertainty. In particular, (i) the impact of UV 
modes k m < k < 0(1 /i/Mpc) for which the perturbative expansion becomes unreliable, and 
(ii) the impact of even smaller scales k > 0(1 h/ Mpc) for which the ideal fluid description 
breaks down. 

Our estimate of effective viscosity and pressure terms in sect. 3.1 is based on one-loop 
expressions and in particular, (rematch contains an integral over the linear power spectrum 
that extends into the deep UV. It is at this point not clear what corrections a s and a v 
receive from these modes beyond the one-loop approximation but one might expect terms 
that resemble the corrections to the power spectrum at low k. The linear density-density 
spectrum P L (q) and its one- and two-loop corrections are depicted in Fig. 2 for z = fj = 0 
when non-linear corrections are largest. It is apparent that the loop corrections are small 
for low k , while they become comparable to the linear spectrum in the region ~ 0.4-1 
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h/ Mpc, say, and finally dominate for even larger wavenumbers.' The analytical estimate 
of the one-loop correction to the low -k spectrum, given by eq. (3.4), involves an integral 
that includes the deep UV region. Nevertheless, the integral of eq. (3.4) is dominated by 
modes k <C 0(1 h/ Mpc) and the contribution from the problematic UV region is small. 
The explanation can be traced to the form of the linear spectrum, which falls off in the 
UV as ~ k~ 3 (up to logarithmic corrections), in combination with the decoupling property 
dictated by total momentum conservation, which ensures that the impact of UV modes 
q on the power spectrum at low k is suppressed by a relative factor ( k/q) 2 [4, 53, 54], 
A similar observation can also be made for the two-loop contribution and in principle for 
higher loops. Nevertheless, the range of wavenumbers which give the dominant contribution 
shifts to higher and higher values with increasing loop order (roughly logarithmically), such 
that they become increasingly sensitive to the UV part of the spectrum. It has been argued 
in [7] that the perturbative series of the power spectrum in the low-A: limit displays the 
behavior of an asymptotic series. A resummation based on Pade approximants suggests 
that the non-linear corrections to the power spectrum at low k are dominated by the one- 
loop result, as one may have naively expected, and correspondingly only mildly dependent 
on deep UV modes. 

It must be emphasized at this point that there is independent evidence from iV-body 
simulations that the effect of UV modes on the low-A: ones through the non-linear mode¬ 
mode coupling is weak [55, 56]. This finding is also consistent with an analytical argument 
that states that small scale structures of size 1/q which are virialized have a milder impact 
on larger scales 1/A: than seen in perturbation theory. In fact, they are suppressed at least 
as (k/q) 3 [54] (see [57] for a recent discussion). Moreover, it was found in [19] that the 
dependence of the power spectrum extracted from coarse-grained IV-body data is consistent 
with the dependence on k m obtained from eq. (3.5) for low k m . 

The impact of very short scale perturbations, for which the ideal fluid picture breaks 
down, on the power spectrum at low k has been investigated in [58]. This was done by 
extracting the second moment of the velocity distribution from IV-body data. The diver¬ 
gence of this velocity dispersion contributes on the right-hand side of the Euler equation, 
and could be interpreted in terms of effective pressure and viscosity as well. Ref. [58] finds 
that the impact on the power spectrum is small, below % level for k < 0.2 h /Mpc at z = 0, 
which would lead correspondingly to small values a s ,a u <C 1 when interpreted as effective 
pressure and viscosity (assuming k m well below 1/i/Mpc). See also [24, 59, 60] for related 
discussions. 

' The spectrum P L appears in the analytical estimate of the one-loop result for the directly observable 
matrix of density-density, density-velocity and velocity-velocity correlations [52] 

( 61 25 \ 

1 | J *?<r$P L (k). (3-1-2) 

This expression is valid at all times, indicating that the one-loop correction (and a fortiori our definition 
of effective viscosity and pressure) diminishes for growing redshift. A similar expression, involving an addi¬ 
tional power of P L (k), holds for the two-loop correction [9, 52]. Our discussion in the previous subsection, 
which was based on the effective propagator, is equivalent to the analysis of the perturbative spectrum 
when projected on the growing mode. 
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Altogether these findings suggest that, for a coarse-graining scale k m below 0(1 h/ Mpc), 
the estimates in the previous section capture the dominant contribution to the effective 
pressure and viscosity, with the largest contribution coming from scales close to k rn . Indeed, 
the integral in eq. (3.5), on which the estimate is based, is dominated by the region near 
the lower limit k m . It is apparent then that k m cannot be taken arbitrarily large, as this 
would move us into a region in which the perturbative expansion fails, making the deduced 
values of viscosity and sound velocity unreliable. At the practical level, we expect that the 
matching scale k m must be taken in the region 0.4 — 0.8 h/ Mpc, with the upper value lying 
at the edge of applicability, and the lower value chosen in a way such that the BAO range 
is well covered. A crucial test of the validity of our assumption is the independence of our 
final results on the specific value of k m within the range discussed above. We shall verify 
this independence in the next section. 

We point out that our construction does not constitute a mere rewriting of standard 
perturbation theory. The crucial difference lies exactly in the way we deal with the UV 
modes. We isolate the effect on the low-A; modes arising from modes near and slightly 
above 1 h/ Mpc, and translate it into effective low-energy parameters (viscosity and sound 
velocity). We expect that these give the dominant effect, while the influence of the deep 
UV is negligible. Nevertheless, it is important to keep in mind that eventually one expects 
small corrections to the effective pressure and viscosity. We then use perturbation theory 
within the low -k theory, cutting off all integrations with an upper limit equal to k m and 
fully taking into account the wavenumber and time dependence of the propagator that 
arises from the effective pressure and viscosity. We expect a quick convergence of the 
perturbative series, a property that we shall verify in the following section. On the other 
hand, standard perturbation theory becomes increasingly sensitive to the deep UV at higher 
loops, a feature that prohibits its convergence. 

4 Results 

Computing non-linear corrections to the power spectrum within viscous fluid dynamics 
is complicated by the wavenumber- and time-dependence of the additional pressure and 
viscosity terms (2.56) in the linear propagation, as well as additional vertices (2.43) and 
(2.45). We use two different methods, which yield results that are in agreement. One of 
the methods turns out to be efficient enough for computing the power spectrum up to 
two loops within the viscous theory, with time- and wavenumber-dependent propagation. 
Before discussing our main numerical results, we briefly describe the methods to compute 
power spectra within the viscous theory, as well as the sensitivity to various assumptions. 
Our numerical results are based on a ACDM model with parameters a$ = 0.79, f i m = 0.26, 
n b = 0.044, h = 0.72, n s = 0.96. 

4.1 Numerical solution 

The first method we use is based on the time-flow approach developed in [14]. It is based 
on coupled evolution equations for the power spectrum and the bispectrum. They follow 
in a straightforward way by taking time-derivatives of the power- and bispectra of density 
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and velocity perturbations, and using the equation of motion (2.55). The non-linear term 
leads to a coupled hierarchy of evolution equations. The time-derivative of the power 
spectrum depends on the bispectrum, and the time-evolution of the bispectrum depends on 
the trispectrum, i.e. the four-point function. Within the time-flow approach, the system 
of equations is closed by neglecting the connected part of the four-point function. The 
bispectrum is then sourced by the disconnected parts, which are quadratic in the power 
spectrum. 

We adapt the method described in [14] in several ways: first, we take the wavelength- 
dependence in the linear propagation into account in the evolution of the bispectrum, as 
described in [57]. Second, we include the additional non-linear vertices (2.43) and (2.45). 
Third, we implement a non-zero initial condition for the bispectrum as described in [57, 61]. 
This is important in order to obtain reliable results when initializing the evolution at a finite 
redshift, as is mandatory in practice. Fourth, we solve the time-flow equations using the 
linear solution for the power spectrum in order to compute the source term in the evolution 
equation for the bispectrum. This is equivalent to a one-loop computation, which however 
takes the time- and wavenumber-dependence of the propagation into account. As was 
shown in [61], this approximation differs only marginally from a self-consistent solution of 
the flow equations for ACDM cosmologies. 

To go beyond one-loop accuracy within the time-flow approach, it is necessary to go 
one order higher in the hierarchy, i.e. to include the (connected) trispectrum. Indeed, 
it was recently shown that including the trispectrum is important to ensure that the re¬ 
sulting higher-order corrections to the bispectrum respect the so-called consistency rela¬ 
tions, which are non-perturbative relations associated with Galilean invariance and with 
the response to local curvature [62]. Unfortunately, solving the time-flow equations with a 
non-zero trispectrum becomes rather cumbersome, even for Einstein-de Sitter and a pres¬ 
sureless, ideal fluid description for which significant simplifications are available due to the 
wavelength-independent propagation [63]. 

Since we are interested in a situation with wavelength-dependent propagation, we opt 
for another approach, which turns out to be more efficient in this case. This allows us to go 
to the two-loop level, taking the wavenumber- and time-dependence in (2.56) into account. 
This second method is a more direct generalization of the standard perturbative approach. 
We expand the density contrast 5 = <£>i and rescaled velocity divergence —0/{Hf) = in 
powers of the initial perturbations at fj = rjo, 

4>a( k, 77 ) = d 3 q 1 ---d 3 q n ( 2vr) 3 5 (3 ) (k- ^ qi)-P’n,o(qi, ■ ■ ■, Qn, ^)<^ qi (%) • • • <W(%) • 

n ^ i 

( 4 . 1 ) 

Inserting this series into the equation of motion (2.55) gives an evolution equation for the 
kernels F n , a , 

n 

(dfjfiab fj))F ^(Ql, • • • ? Qn? V) ~ ^ ^ '7abc(Q 1 “1“ ’ ’ ’ Qm? Qm+1 ’ ’ ’ H - Qn) (4.2) 

m= 1 

X 1 > • • • ) Qrrtj d)Fn —m,c (*lm+l > • • • ) Qru fj) ) 
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where the right-hand side is understood to be symmetrized w.r.t. arbitrary permutations 
of the q 2 , and k = ]>Tqj. When neglecting the pressure and viscosity terms, the solu¬ 
tion coincides with the standard Einstein-de Sitter kernels [5], F n> i| a „=a s =o = e nv F n and 
Fn, 2 \a„=a s =o = e nf ^G n . For non-zero pressure and viscosity, the time-dependence does not 
factorize and it is hard to do the time-integration analytically. We therefore solve the dif¬ 
ferential equations numerically. Since the effective pressure and viscosity become negligibly 
small at early times, we initialize the kernels with their Einstein-de Sitter values at an early 
initial time chosen to be f) = —4 (we checked that our results are stable when varying the 
starting time). The equation (4.2) allows for a recursive determination, starting from the 
linear solution (n = 1). This solution is then used to determine the kernels with n = 2, etc. 
To obtain the power spectrum at one-loop, one needs to go up to n = 3, and at two-loop 
up to n = 5. 

The power spectrum is then obtained by the same expressions as in standard pertur¬ 
bation theory (i.e. Pi-i oop = P 22 + 2 Pi 3 , P 2 -I 00 P = P 33 + 2 P 15 + 2P 24: ) [5, 52], except 
that the kernels are given by the time-dependent solutions of (4.2) instead of the Einstein- 
de Sitter kernels. For an efficient implementation it is essential to compute the kernels 
only for the wavevector configurations that are actually required for the loop integrals, 
e.g. P 5 ja (k, q, —q, p, —p, fj) for the two-loop contribution Pi 5 , along with the lower-order 
kernels required for the recursive solution. Using the manipulations of the integrand and 
the algorithm described in [7], this method allows us to compute the power spectrum up 
to two loops within the viscous theory at an affordable numerical cost (of order a couple 
of minutes per external wavenumber k on a single desktop machine). We checked that the 
one-loop results obtained with both methods are in excellent agreement. 

4.2 Pressure versus viscosity and extra vertices 

The matching procedure described in subsection 3.1 allows us to obtain separate expressions 
for the parameters a v and a s characterizing the effective viscosity and pressure, respec¬ 
tively. However, only their sum enters the growing mode of the linear solution in the low-k 
limit. Nevertheless, the non-linear contributions are sensitive also to decaying modes inside 
the loops, and to wavenumbers for which the low-k limit is not applicable. Therefore, one 
may wonder how sensitive the power spectrum is to changing the ratio a u /a s while keeping 
their sum fixed according to (3.10). In fig. 3 we show the dependence on this ratio for a 
matching scale k m = 0.6 h/ Mpc at z = 0. While the linear solutions differ for k > 1 h/ Mpc, 
we find that both the linear and the one-loop result for the power spectrum agree to better 
than 1% for k < k m . Thus, while pressure and viscosity lead to qualitatively different 
behaviour at large k, their effect is degenerate within the weakly non-linear regime. This 
is consistent with the analytical expression for the linear propagator in eq. (2.39) and in 
a similar form, this feature has also been observed within a perturbative treatment of 
effective pressure and viscosity terms in [ 22 ], 

Apart from a more complicated propagator, viscosity and pressure in general also 
induce additional non-linear interaction vertices (2.43) and (2.45). In fig. 4 we show the 
one-loop solution obtained when taking these vertices into account, compared with the case 
where only the standard vertices are kept, again for k rn = 0.6/i/Mpc and at z = 0. We 
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P d5 ( k,z=0), k m =0.6h/Mpc 



Figure 3. Power spectra in linear and one-loop approximation for a viscous fluid. The lines 
(dots) show the linear (one-loop) solution obtained for various ratios of viscosity and pressure 
terms. Their sum is kept constant. The differences between the various results are below the % 
level for k < 0.6/i/Mpc. For comparison, the black solid line shows the linear spectrum obtained 
in the pressureless, ideal fluid case. 

P S6 ( k,z=0), k m =0.6h/Mpc 



Figure 4. One-loop power spectrum for a viscous fluid at z = 0 for two different approximations 
of the non-linear coupling terms. The blue dots correspond to the result when taking the standard 
vertices (2.42) and (2.44) and the extra vertices (2.43) and (2.45) into account. The green squares 
show the result when neglecting the extra vertices. 
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find that the difference is below the 1%-level for k < k m . We therefore neglect these extra 
vertices in the following. This finding also justifies why it is not crucial to resolve possible 
ambiguities in the precise magnitude of these extra vertices, as discussed in subsection 2.3. 


P 5S ( k,z=0), k m =0.6h/Mpc 


P 6S ( k,z=0.375), k m =0.6h/Mpc 




P w (k,z=0.833), k m =0.6h/Mpc 


P«(k,z=1.75), k m =0.6h/Mpc 




Figure 5. Power spectrum obtained in the viscous theory for redshifts z = 0, z = 0.375, z = 0.833 
and 2 = 1.75, respectively, normalized to the conventional linear spectrum Pu n ,ideal (i.e. assuming 
an ideal fluid). The open (filled) circles show the one-loop (two-loop) result computed in the viscous 
theory. The solid blue line is the linear spectrum computed in the viscous theory, and the red points 
show results of an IV-body simulation [64], 


4.3 Comparison with IV-body simulation and dependence on the matching 
scale 

In fig. 5 we show the matter power spectrum at various redshifts, computed for a viscous 
fluid. All curves are normalized to the linear spectrum obtained for an ideal, pressureless 
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fluid. The solid lines show the linear spectrum in the viscous theory, which is suppressed 
compared to the ideal fluid for large k and small z, as expected. The open (filled) circles 
show the one-loop (two-loop) result, taking the time- and wavelength-dependent viscosity 
and pressure terms into account, and including modes with q < k m in the loop integrals. 
The viscosity and pressure are determined from the matching procedure described in sub- 
section3.1 with k m = 0.6/r/Mpc. 
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Figure 6. Dependence of the power spectrum on the matching scale k rn , for z = 0 and z = 
0.375. The open (filled) symbols show the one- (two-) loop results as in fig. 5, but for various 
values of k m . The lines show the linear results. The numerical values of the effective pressure 
and viscosity coefficients obtained from the matching condition (3.10) are given by a u + a s = 
(1.02,1.34,1.60,1.80) • Dl(z ) 2 for k m = (0.4,0.6, 0.8,1.0) /i/Mpc, respectively. 

For comparison, the red points show corresponding results from a large-scale IV-body 
simulation [64], The two-loop results agree reasonably well with the IV-body data for 
k < 0.2 — 0.25/i/Mpc at low redshift, and the agreement improves at higher redshift. It 
must be emphasized that no parameter adjustment has been made in order to reproduce 
the IV-body results. 

Since the full result should be independent of the choice of k m , any residual dependence 
can be taken as a measure of the theoretical uncertainty. In fig. 6 we show the dependence 
on k m in the range 0.4 — 1 h/Mpc. The variation of the two-loop results is small within 
the range for which the results agree well with IV-body data. The size of the two-loop 
correction increases with k m , i.e. also the perturbative uncertainty becomes larger. For 
the largest value k m = 1 h/Mpc, the deviations start to become sizeable, as expected from 
the previous discussion. As mentioned above, the residual dependence on the matching 
scale can be taken as an estimate for the theoretical uncertainty. Considering a variation 
over a factor of about two, in analogy to common practice within quantum field theoretical 
computations, indicates that our results are robust at the % level for k < 0.2 h/Mpc 
at z = 0. As discussed in section 3.2, the remaining dependence can be attributed to 
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uncertainties in the determination of the effective pressure and viscosity. In addition, 
deviations with respect to IV-body simultations are expected to occur for k m > 1 h /Mpc 
or k > 0.2/i/Mpc because either the coarse-graining starts to become sensitive to scales 
for which the usual fluid description becomes inaccurate, or higher-order non-linearities 
become important, respectively. 


P w (k,z=0)/P w _body, viscous theory P«5a(k,z=0)/P w _body. SPT with cutoff 




Figure 7. Comparison of results for the power spectrum obtained within the viscous theory 
(left) and in standard perturbation theory (right) at one- and two-loop level. In this figure we 
normalized the power spectra to the IV-body result [64] at z = 0. The grey band corresponds 
to a ±1% deviation, plus an estimate for the error of the IV-body results based on the boxsize 
L = 7200Mpc /h and number of particles N = 6000 3 . The left plot shows our one- and two-loop 
results within the viscous theory for various choices of the matching scale fc m , as in fig. 6. For 
comparison, the right plot shows one- and two-loop results in standard perturbation theory (shown 
as dashed and dotted lines, respectively), computed with various values of an ad-hoc cutoff A 
(coloured lines), as well as in the limit A —> oo (black lines). 

Importantly, also the size of the two-loop correction relative to the one-loop is smaller 
by a factor of a few as compared to standard perturbation theory. This is partially due 
to cutting off loop integrations at k m , but can also be attributed to the reduced UV 
sensitivity within the viscous theory. This can be observed in fig. 7, where we compare our 
results (left plot) in the viscous theory with standard perturbation theory results (right 
plot). For the latter, we use an ad-hoc UV cutoff A chosen to coincide with the matching 
scale k m . It is apparent from fig. 7 that the two-loop result within standard perturbation 
theory depends strongly on the cutoff A. In contrast, as discussed earlier, the one- and 
two-loop results within the viscous theory are significantly less sensitive to the matching 
scale. We are aware of two effects that can contribute to a reduced UV sensitivity: first, 
as discussed in subsection 3.2, the effective pressure and viscosity are dominated by modes 
close to k m , i.e. the impact of very short scales on modes with k < k m is actually small (in 
line with the findings of [54-56], as discussed before). Second, the effect of pressure and 
viscosity itself leads to a damping of UV modes, as discussed in detail in subsection 2.2 
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for the linear solution. This damping effect propagates also to the non-linear solution, in 
particular through the wavenumber-dependent terms in (2.56) appearing on the left-hand 
side of (4.2). For the second mechanism to be effective it is necessary to compute the 
kernels within the viscous theory, rather than expand them in the pressure and viscosity 
contributions. The reason is that the latter strategy would lead to contributions that grow 
with powers of the wavenumber, which artificially enhances the UV sensitivity, instead 
of reducing it. On the other hand, the numerical treatment pursued here captures the 
damping effect which sets in for modes close to k m . 

We conclude that the two-loop results are robust with respect to a variation of the 
matching scale over the range ~ 0.4 — 0.8 h/ Mpc, and display an improved agreement with 
iV-body data as compared to the two-loop computed in standard perturbation theory. 
In addition, our results suggest that the convergence of the loop expansion is improved 
compared to the pressureless, ideal fluid case. This is an important property because it 
reduces the perturbative uncertainty. 

5 Conclusions 

We explored here a description of dark matter as a viscous fluid with non-vanishing sound 
velocity. For scalar perturbations in the regime of dark matter domination one finds that 
viscosity has the tendency to slow down gravitational collapse and structure formation in 
the linear regime for modes with large wavevector. However, in contrast to the effect of 
pressure, viscosity does not stop structure formation completely, and it does not wash out 
structures. 

These properties are independent of whether viscosity and pressure have microscopic 
origin in terms of a non-trivial velocity distribution of the constituents of the dark matter 
fluid, or whether these terms arise on an effective level from the non-linear behavior of 
short wavelength modes. It is this latter possibility that we explored in more detail. 

More specifically, we obtained concrete expressions for the effective viscosity and sound 
velocity coming from short modes still within the fluid approximation by matching the lin¬ 
ear propagator of the viscous fluid to results for the propagator beyond linear approximation 
in standard perturbation theory. The resulting values for viscosity and sound velocity de¬ 
pend on the cosmological time and on the spectrum of inhomogeneities. Physically, these 
terms account for the leading effect that modes with wavelengths below a coarse-graining 
scale 2ir/k m have on modes above that scale. 

For the modes with wavevectors k < k m we numerically solved the viscous fluid equa¬ 
tions in linear, one-loop and two-loop approximation. The resulting density-density spec¬ 
trum in the BAO range up to k < 0.2h/Mpc agrees on the percent level with the results 
from ./V-body simulations. This level of precision is consistent with the expected size of 
corrections from deep UV modes for which the fluid approximation becomes inaccurate. 
The theoretical error can be estimated by the sensitivity to the matching scale k m . We 
find that the two-loop results within the viscous theory are robust at the percent level to a 
variation within the range 0.4 — 1 h/ Mpc, in marked contrast to the cutoff dependence of 
the two-loop result within standard perturbation theory. Moreover, the effective pressure 
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and viscosity lead to improved convergence of loop integrals in the ultraviolet so that one 
may hope for an improved convergence of perturbation theory also at higher orders. 

In this sense, the theory with effective viscosity and pressure offers a convenient calcu- 
lational framework for understanding inhomogeneities in the cosmological fluid at interme¬ 
diate wavelengths. This set-up may be rather efficient in order to include also the leading 
effect of other phenomena that become relevant at short wavelength such as rotational 
(vector) degrees of freedom or heat transport due to interactions of the baryonic matter 
component. Finally, it would be intriguing to study the viscosity and pressure arising from 
kinetic motion and interactions in the dark matter sector in order to better understand its 
microscopic properties. 
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